Waves on Icicles 
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Icicles with wave patterns on their surfaces can sometimes be seen hanging from roofs of buildings. 
Surprisingly, most of these wave patterns are at intervals of about 1 cm. The reason for this 
uniformity of interval has not been clarified. Here we show a formula to explain this remarkable 
phenomenon by introducing a new instability theory. This theory is given by thermal diffusion in 
thin water layer streams flowing along the icicles. The streams change the temperature distribution 
and control waves of short wavelengths. The specific wavelength (about 1 cm) can be determined 
by Laplace instability of the heat field in the atmosphere and by the thermal diffusion effect in 
thin-layer streams. 

PACS numbers: PACS number(s): 47.20.Hw, 81.30.Fb 
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Icicles are often seen hanging from roofs of buildings 
in winter season in snowy regions. They are almost 
columnar in shape and their radii (or thickness) become 
smaller in the gravitational direction. The precise inner 
morphology of icicle is much more complicated and in- 
teresting like having hollows, which is shown by Maeno 
et.al. jTj. But in this article we discuss only the surface 
instability of icicle as the periodically modulated thick- 
ness of an icicle along its longitudinal direction. The 
intervals of wave (wavelength) appearing on icicles are 
almost uniform (about 1 cm), and do not depend on the 
length or thickness of the icicle or on the external tem- 
perature. Furthermore, such waves appear only on the 
surface of icicles that are enveloped by a layer of water 
flow of about 0.1 mm in thickness. The surface of a grow- 
ing icicle changes from flat to wavy in a few hours [^|. 
Formation of such waves is also discussed qualitatively 
by Maeno et.al. |j] as the conjecture. But its scenario 
can not explain both the value of wavelength and its uni- 
versality. The mechanism to construct the wave here we 
show, is different from above one. 

This kind of phenomenon is usually called surface in- 
stability. In this short article, we will discuss both the 
instability of the surface of an icicle and factors that de- 
termine the wavelength appearing on the surface of an 
icicle. In 1997 Matsuda has given the experiment to con- 
struct the wavy ice in which the physical situation is 
the same as icicle. He has made a gutter inclined at 
8 degrees and kept the static water flow on it. He put 
the apparatus in cold room with temperature -8 degrees 
centigrade. Then the wavy ice has appeared under the 
water flow and its wave length is about 8 mm at = 7r/2, 
which agrees with the case of icicle Q . This means the 
waves appearing on icicle is induced by the interaction 
between thin water flow and thermal diffusion. 

In the study of crystal growth, surface instability is 
usually explained by Mullins-Sekerka's theory ||, which 
is based on two observations: Laplace instability and 
Gibbs-Thomson effect. However, the mechanism un- 
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derlying the instability that occurs on icicles cannot be 
explained by this theory. The reason for this is that 
the Gibbs-Thomson effect only applies to wavelength 
of about 10 micrometers; i.e., wavelength much shorter 
than those on icicles (about 1 cm). Furthermore, Laplace 
instability does not occur in a thin liquid layer; it only 
occurs when the thickness of the layer is wider than the 
wavelength. Thus, Laplace instability could only occur 
on icicles if the thickness of liquid layer were larger than 
the wavelength of 1 cm, which is never the case. 

Let us then consider the characteristics of a fluid. 
There are various instability theories in fluid mechan- 
ics. Benney's liquid film has often been used to study 
the characteristics of a thin water layer . This theory 
explains why waves appear on the surface of thin water 
flow on a sloping road on a rainy day. This phenomenon 
is caused by gravity, viscosity and surface tension. How- 
ever, these are not static waves but moving waves. In 
the case of icicles, waves move down an icicle two-times 
faster than the average velocity of a fluid. Thus, the for- 
mation of waves on an icicle can not be explained only 
on the basis of characteristics of such hydrodynamical 
waves. 

We therefore consider not only the fluid characteristics 
which flows along an icicle but also the thermal diffusion 
effect to express the static wave on the surface of ici- 
cle. Our program is the following. First, we solve the 
hydrodynamics of thin-layer fluid, which is not used in 
the usual Mullins-Sekerka theory. Next, we consider the 
thermal diffusion in thin-layer fluid and in the air sur- 
rounding an icicle. The thermal diffusion in air is very 
important since Laplace instability cannot occur in a wa- 
ter layer. As a boundary condition, the heat flow should 
be conserved on the boundary between air and the wa- 
ter. The thermal properties in water fluid and in air give 
all the necessary information on the formation of surface 
of icicles. 

Let us start with the case of an icicle having a bound- 
ary surface (ice-water) with a small perturbation. Wc 
add a perturbation to the radius R of an infinitesimally 
long columnar ice (model of icicle), 
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R-> R + S(t) sin kx, 



(1) 



(x-axis is along the length of the icicle) and we consider 
the time development of S(t). We define the y-axis as 
perpendicular to the ice surface, and use the variable 
y = r — R with | y \<< R. The thin water layer on 
this surface has mean thickness ho ~ 0.1mm. For con- 
sidering such a film fluid, it will be easier to have a fluid 
dynamical model as thin water layer on a inclined gutter 
just like Matsuda's experiment [g). (x-axis is taken along 
slope, and y-axis is taken normal to slope surface.) This 
simplification is well established since ho <C R- This sit- 
uation for fluid mechanics is usually discussed for flat 
gutter, and shown the specific surface velocity 
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with parabolic velocity distribution, 

«* = CT(2f -(f) 2 ), 
ho h 



(2) 
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where v is the viscosity (1.8 x 10~ 6 m 2 /s) and g is the 
gravitational acceleration. (For this calculation, see for 
e.g. the exercise in Landau's text book [Q.) If we observe 
mean surface velocity U, we obtain ho and flow quantity 
Q = 2irRUho, where U = j~ J Q ha v x dy = 2U/3 is the 
mean velocity. In Matsuda's experiment || , he used the 
flow Q = 160ml /hr and width of gutter W = 3cm (W 
corresponds to 2irR), where the wave appears most clear. 
This gives Uh = 1.48 x 10~ 6 m 2 /s. From U = 2U/3 

and U = ^f, we get U = 2.4 x 10~ 2 m/s with h = 
0.93 x 10 _4 m. On the other hand, his direct observation 
of surface mean velocity by using test fine grain gives 
U = 4 X 10~ 2 m/s. So we use U = (2.4 ~ 4) x 10" 2 m/s 
with ho = (0.93 ~ 1.21) x 10~ 4 m as the experimental 
surface velocity and thickness of water layer. 

Our problem with small modulation S sin kx on bound- 
ary can be expressed as perturbation to above solution. 
To do this we make nondimensional coordinate a;* and 

as 

a* = kx, y* = y/h , 

and we put off * in a while. In this notation solid surface 
is given as y — 77 sin a;, where r\ — 5 /ho. The fluid in such 
a thin-layer has boundary conditions as follows. The 
fluid velocity equals to zero at the ice-water boundary 
due to viscosity, and the gradient of fluid velocity equals 
to zero on the air-water surface due to the vanishing of 
shearing stress. 

By solving static Navier Stokes equation |5| with above 
mentioned boundary conditions and with long wave- 
length approximation ji = kho ~ 10~ 2 , we obtain the 
hight of air-water boundary surface in nondimensional 
form given as, 



which means the thickness of water layer is uniform at 
lowest order of long wavelength approximation. (In first 
order approximation, we have modulation in its thick- 
ness but it is not thinner at protrusions. The thickness 
modulation is like cos kx against the icicle modulation: 
sinfca;. But these effects are irrelevant. ||) The stream 
function of water is obtained as 

V^* = "^(^ ~ ^sinx) 3 + (y - r/smx) 2 , (5) 
where ifi* is also the non-dimensional stream function as, 
tp = hoUtp*, 

and we put off * again. The form of such stream function 
is not surprising, and it means just the parabolic form of 
velocity fields, where the velocity vanishes at modulated 
solid-water boundary. 

The next task is to solve the thermal diffusion equation 
with a flow background as, 
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where D = — ~ 1.3 X 10 7 m 2 /s is the thermal diffusion 
constant and we used non-compressive condition. We 
rewrite above form by using stream function tj;, neglect- 
ing higher order in parameter fi = kho as long wave- 
length approximation, and under the static thermal dis- 
tribution condition. (All variables are written in dimen- 
sionless form except temperature.) 



d 2 T 
dy 2 



dT dip dT dtp 
dx dy dy dx 



(J) 



Note that the dimensionless parameter a = ^XJho/D 
is proportional to wave number k. The value of a is 
about 2 experimentally. By using boundary condition 
on solid- water surface (y = rysina;)as 
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Q(SW) = -k^7 \ v=V sinx= -m(x), (9) 



T(SW) = T M = const, 
dT 
dy 



we obtain the quantities on air-water surface 
(y = 1 + vsinx) by solving differential equation. 
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D 2 ]a{x), (10) 
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D 2 ]a{x), (11) 



nsmx, 
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where Tm is the melting point, Q's are heat flows, a(x) 
is an unknown gradient of temperature on solid surface, 
and D means ad/dx.The higher order terms of a are 
negligible since their coefficients are effectively small. 

It should be noted that the temperature of the air- 
water surface is different from the one in the environ- 
ment. The temperature in atmosphere changes greatly 
near the surface of an icicle. Thus, we must also consider 
the thermal diffusion in air. Since we consider the static 



3 



diffusion, our equation becomes Laplace equation. We 
write it in cylindrical coordinate as 

d 2 Id Id 2 d 2 

[^ + -^ + ^+^]nr,M)=o. (12) 

Note that our previous coordinate y is related to r as 
r = R + y. 

We suppose the existence of axial symmetry, and we take 
dT /dO = 0. Therefore we work with 
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(13) 



Since we have oscillating surface in x-direction, it is nat- 
ural to construct the solution in the form, 



T(r, x) = f(r) + g(r) sin(fcx + 

where f(r) satisfies 

^d 2 + 1 d q 
dr 2 r dr 

and g(r) satisfies 

dr z r dr 



(14) 



(15) 



(16) 



Then the solution is given by using dimensional vari- 
ables x,y = y — 1~iq,R = R + 1iq as 

T(x,y) = A + B\og(l + y/R) 

+ CK Q (k(R + y))sm(kx + 4>), (17) 

where A,B, and C are the coefficients and Ko is the modi- 
fied Bessel function. Last term is induced from boundary 
condition at water surface. The reader might think the 
appearance of logarithmic term strange, since it diverges 
at large r. But the appearance of such a term is natu- 
ral for infinitesimally long axially symmetric source. In 
our system the length of icicle is finite. Therefore this 
solution is valid just around the icicle, and far from ici- 
cle source we must treat the source just like point and 
we need to connect two solutions at r ~ L (L is the or- 
der of length of icicle). In this way the coefficients A,B 
and C will be determined. But we have another physi- 
cal observable, such as, mean growth rate of icicle which 
determine coefficient B. The information of temperature 
at infinity is included into this observable. 

At y = 8 sin kx (AW surface) , temperature and ther- 
mal flow are given by above solution by using taylor ex- 
pansion. 
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A + [CK Q (kR) cos (j) 8] sin kx 
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+CKo(kR) sincficoskx, 
LV - [K CkK' (kR) cos0- 
—KoCkK' (kR) sin</>cos kx, 



(18) 
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V is the mean growth rate of icicle defined by 

V =<Q> /L, 

with L the latent heat, and < > means the average 
along x-direction. This value determines coefficient B 
as, 



B = — 
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kq is the thermal conductivity of air. We treat C as 
order S and we took up to first order in 5. 

Note that the thermal-diffusion inside an icicle is not 
considered. Since the surface of icicle is the same temper- 
ature as melting point, the temperature inside an icicle 
can be considered as uniform. We now give continuity 
condition of temperature and heat flow conservation on 
the water-air boundary surface. By replacing (1C) and 
(11) in dimensional form, and compare with (18) and 
(19), all the coefficients are determined as well as the 
function a(x). The growth rate of solid is determined 
from this function as 

v(x) — —na(x)/L = V + bs'mkx + ccoskx, (20) 



where b and c are obtained some coefficients. We can ob- 
tain the amplification rate of fluctuation by 5/8 = b/8. 
We have two mixed effects for the growth rate of fluctua- 
tion. One is the diffusion property of temperature in the 
atmosphere. The other is the effect of fluid for thermal 
diffusion in thin-layer. In total, we obtain the amplifica- 
tion rate of fluctuation under the condition kR >> 1/2. 
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where a cx k , and V is the mean growth rate. This 
equation is interpreted as 

R.H.S = (Laplace instability) 

x (fluid effect for thermal diffusion). 

Therefore, thermal diffusion in atmosphere enhances 
fluctuation of shorter wavelength, but fluid effect for 
thermal diffusion tries to suppress such fluctuation of 
shorter wavelength. From these two effects, amplifica- 
tion rate has a peak value. The maximum amplification 
factor is given at definite wave number. 

O^rnax ~ 2.2. 

The corresponding wavelength is obtained by using the 
definition of a as, 

~ Da 

The experimental data: 

U = (2.4 - 4) x l(T 2 m/s, h = (0.93 - 1.21) x 10~ 4 m 
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with a max ~ 2.2. gives A = 5mm ~ 13mm, which well 
agrees with observed wave length 8mm. 

Let us discuss our obtained result by comparing with 
the discussion given by Maeno.et.al.Q. Our obtained re- 
sult comes from two important effects. One is the ther- 
mal diffusion in air and its Laplace instability another 
one is the effect of fluid for thermal diffusion in thin 
water. The fluid effect shuffles the static temperature 
distribution and so the higher frequency wave is sup- 
pressed. On the other hand Maeno et.al. jjj have given 
the conjecture on the formation of ribs on icicle. This is 
from three points as follows. 

1: If we had an protrusion accidentally the water flow 
on it may be thinner and cooled much faster. 2: The 
protrusions are more exposed in the surrounding cold 
air and giving rapid growth of ribs. 3: Since the water 
flow have larger viscosity the water tends to be stagnant 
just upstream of ribs and giving slower growth of icicle. 

The first and third one are irrelevant in long wave- 
length approximation which we used, and the thickness 
of water layer is almost uniform in our explicit fluid me- 
chanical calculation. This will be well understood by 
the difference between thickness 0.1mm and wavelength 
10mm. The second one may correspond to the Laplace 
instability in air which we used. But this is not enough 
to express the wave properties. Really all of above three 
remarks enhance shorter wavelength fluctuations, and 
do not determine specific wavelength. So we have used 
not only the property of water thickness, but also the 
velocity distribution of water flow which affects thermal 
diffusion. This effect breaks the Laplace instability at 
shorter wavelength, and we obtain peak growth rate at 
definite wavelength as the quantitatively good result. 



The wavelength depends on the fluid flow 
Q = 2TrRUho on icicle surface: 

A roax = ^ = J^(-) 1/a (§r /3 , (22) 

L>OL m ax t)a max giT ZH 

where R is the mean radius of the icicle. The tempera- 
ture of the environment changes the mean growth rate of 
solidification, but wavelength is independent of it. The 
wavelength depends on the ratio of Q to R but not only 
on Q. We then have a universal wavelength by assuming 
that Q is proportional to R. The fluid effect that sup- 
presses short wavelength fluctuation is somewhat simi- 
lar to the Gibbs-Thomson effect. But their origins are 
quite different, and they work in different length scales. 
Further experimental test will be necessaryto check the 
conclude validity of this theory by using But here 

we have shown how the wave can be presented on icicle 
as instability theory. Our proposed instability theory is 
quite new and may be applied to a wide range of phenom- 
ena. For example, the similar phenomenon occurring on 
stalactites can be explained in the same way by changing 
thermal diffusion to solution diffusion. 
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